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ABSTRACT 

The finite element method (FEM) is a numerical algorithm for solving second 
order differential equations. It has been successfully used to solve many problems in 
atomic and molecular physics, including bound state and scattering calculations. To 
illustrate the diversity of the method, we present here details of two applications. First, 
we calculate the non-adiabatic dipole polarizability of H 2 by directly solving the first 
and second order equations of perturbation theory with FEM. In the second application, 
we calculate the scattering amplitude for e-H scattering (without partial wave analysis) by 
reducing the Schrodinger equation to set of integro-differential equations, which are then 
solved with FEM. 


I. INTRODUCTION 

The finite element method (FEM) is a powerful numerical tool for solving 
differential equations, including eigenvalue problems [1]. FEM utilizes a piecewise 
interpolation scheme, in which the unknown function is approximated locally by a simple 
polynomial. Although the method was originally developed to solve problems in 
structural mechanics, Frank Levin was one of the first to recognize that FEM could be 
used to study few-body systems [2], For bound states, FEM has been used to calculate to 
high accuracy the energy and structure of three-body Coulomb systems with arbitrary 
masses [3]. Another important application is the study of atoms and molecules in strong 
external fields [4]; when the wave function is strongly distorted, the piecewise 
interpolation approach is often superior to standard hydrogenic or Gaussian basis set 
expansion. FEM can also be used to study collisions, where the complicated boundary 
conditions associated with scattering can be treated in a straightforward fashion. 

Accurate phase shifts and inelastic cross sections have been calculated for e-H collisions 
for 0 < L < 3 [5]. For the Temkin-Drachman Retirement Symposium, I have selected 
two FEM examples which reflect the research interests of the honorees. 

Richard Drachman has contributed greatly to our understanding of long-range 
interactions. In his series of papers on the Rydberg states of Helium [6] and Lithium [7], 
he has provided a rigorous theoretical description for Rydberg atoms based on an 
effective polarization potential. An important extension of this work has been the 
formulation of an effective polarization potential for the Rydberg electron of H 2 ; 
microwave spectroscopy of high Rydberg states provided a mechanism for determining 
the multipole moments of the ionic core, including the static dipole polarizability [8]. 
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The measurement of a s for motivated several groups (including Drachman and 

Bhatia [9]), to calculate the polarizability of Hj without invoking the Bom- 
Oppenheimer approximation. In section II, we present a description of a FEM 
calculation of the dipole polarizability of Hj . 

Aaron Temkin’s contributions to scattering theory have had a major and lasting 
impact on the field. With his method of polarized orbitals [10], he was the first to include 
the effects of polarization and exchange in the ansatz for the wavefunction. In 1962, he 
introduced the now famous Temkin-Poet model [11], For the past eight years, I have had 
the privilege to collaborate with Aaron, pursuing a new approach to scattering that does 
not use partial wave analysis. The scattering amplitude is calculated directly by solving a 
set of coupled integro-differential equations. Section III summaries our progress to date 
and outlines our plans for the future. Atomic units are used throughout. 


II. NONADIABATIC DIPOLE POLARIZABILITY OF H; 

In the late 90’s, experiments on the Rydberg states of H 2 provided a mechanism 
for determining the static dipole polarizabiltiy of Hj to high precision [8], At that time, 
the only theoretical calculation for a s employed the Bom-Oppenheimer approximation 

[12] , Given the accuracy of the new results, it was not surprising that there was a 
discrepancy between the experimental value a s =3.1681(7) and the theoretical value 

« s bo = 3.1713 on the order of m e / m p . This breakdown of the Bom-Oppenheimer 

approximation motivated several groups to attempt a non-adiabatic calculation of the 
dipole polarizability. 

The static dipole polarizability a s is defined in terms of the second order 
correction to the energy due to the presence of an external electric field s : 

E (2) =-\a s e 2 =('F (1) (F,j?)|(l + S)s r | 'F (0) (F,R)) (1) 

where (1 + S) = . In the variational approach, the first order correction to the wave 

function is expanded in a basis set that includes nuclear and electronic states. Using 
FEM, one can solve directly the first and second order equations of perturbation theory 

[13] . The first step is to carry out the frame transformation that reduces the number of 
variables in the problem. 

Frame transformation 

In the space-fixed laboratory frame ( x' , y' , z' ) , the electric field is aligned with 
the z' -axis (see Fig. 1). After separating out the center-of mass motion, the (field-free) 
Hamiltonian for the relative motion is given by 
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Figure 1. Hj in the space-fixed (laboratory frame). £ indicates the direction of the 
external electric field. 
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and ju = m p / 2 . 

The vector R lies along the intemuclear axis and r is the vector from the center of the 
intemuclear axis to the electron. The Hamiltonian commutes with L 2 and L z , and in 

general, the non-adiabatic wave function v f / (L,M';r,i?) depends on all six coordinates. 

To simplify the problem, we perform a rotation 9?(0' , 0' ,0) which leaves the 
intemuclear axis aligned with the new z-axis. It appears that we have eliminated two 
degrees of freedom, since the wave function is now a function of only r and R . But 
there is a price to pay for this frame transformation. The Hamiltonian does not commute 
with L z ; M is not a good quantum number and the Hamiltonian is not diagonal in the 

basis spanned by the eigenstates of L 2 and L z . The electric field, which appears in the 
matrix element of Eq. (1), is now a function of the Euler angles 

s = £•(- sin ©' x + cos 0' z) . (4) 
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Figure 2. H 2 + in the body-fixed frame, s indicates the direction of the external 

electric field. 


Despite these complications, it is still desirable to work in the body-fixed frame. 
Ultimately, we must reduce the problem to a solution of a set of differential equations in 
three variables if we are to apply the FEM. 

The space-fixed wave function (with ‘good’ quantum numbers L, M’) is a linear 
combination of the body-fixed wave functions, 

'V SF {L,M'-r,R)= (5) 

M=-L 


where D L MM , ,Qi) are the coefficients of the irreducible representation associated 
with the rotation 91 (<!>' , 0' ,0) . The Hamiltonian in the body-fixed frame is given by 
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Perturbation Theory 

We now return to the evaluation of the matrix element in Eq. (1). It is clear that 
we must start in the space-fixed frame, where the ‘good’ quantum numbers associated 
with the unperturbed ground state are L = 0 and M' - 0 ; from the Wigner-Eckart 
theorem, we know that the first order correction to the wave function must be a state with 
L = 1 and M'= 0 . Therefore, we re-express Eq. (1) more precisely as 

E (2) = -}a s £ 2 = (¥$ (1,0; r,R) |(1 + 8)e-r\ ^(0,0 ;r,R)) (8) 

where '¥ ( S ° F ) and V P^ ) are found by solving the zeroth and first order equations of 
perturbation theory: 

H ^ ( 0,0; r , R) = E (0) ^ ( 0,0; r, R) (9a) 

(/f-£ (0) )'P^(l,0;F,^) = -(l + ^)^ r ¥£>(1,0 -,r,R). (9b) 


For the special case L = 0 , M'= 0 , there is a one-to-one correspondence between 
the space-fixed wave function and the body-fixed wave function 

4^>(0,0; r,R) = ^°;(0,0;F,i?) . (10) 

The body-fixed Hamiltonian is diagonal since L ± v Fa° ) (0,0;r,7?)=0. Furthermore, the 
wave function is independent of the electronic azimuthal angle <p , and the resultant 
differential equation (in three variables) 


' 1 

< 

p 2 r +\lI 

+ 

fl 1 1 

p 2 + -jL 2 

\ 2 r 

L‘ R 2 J 


U 8 jur J 

r 


+ V(r, R) - E m \ ^ (0,0; r, 9, R) = 0 (11) 


/( 0 ) 


can be solved with FEM for the ground state energy and wave function. We obtained 
£ (0) =-0. 597 139 055(8). 

In order to obtain the first order correction to the wave function, we need to solve 
Eq. (9b) by transforming it to the body-fixed frame. The zeroth order wave function 
v Fj“ ) is given by Eq. (10) and the dipole interaction is given by 

s r = £ , r(cos0’cos#-sin0'sin0cos^). ^ 
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Using Eq. (5), the first order correction to the wave function V P^ ) is re-expressed in the 
body-fixed frame as 


^ (1,0; F, R) = cos©' 4^ (1,0; F, R) + ^sin 0' [4^ (1,1; F, R) 


/O) 


(1)/ 


■^>(l-l;F,i?)]; 


( 13 ) 


the (body-fixed) Hamiltonian is given by Eq. (6). 

The resultant coupled equations depend on the five variables ©' ,r,6,(/> and 
R . Equating terms that multiply cos 0' and sin 0' , we can eliminate the dependence on 
the Euler angle 0' . We are left with three coupled differential equations for the body- 
fixed wave functions 4^ (1, M;r,R ) , M = 0,±1 . The dependence on the electronic 
azimuthal angle (j) can be obtained analytically; one can show that the wave functions 
must be of the form 

'Ol frr,R) = £f{r,0,R) (14a) 

'CtUl ■r,R) = ±^sg(r,9,R)e ± '* . (14b) 


Thus we have reduced the problem to a set of two coupled equations in three variables: 


[H m -E (0) ]f(r,e,R)+^(^ + cotd)g(r,8,R) = -rco S 8^(^d,R) 

pR \38 J 
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(16) 


Eqs. (15a) and ( 1 5b) are solved with FEM. Once we know f(r,8, (f>) and g{r,8, (fi) we 
can construct the first order correction to the wave function and evaluate the matrix 
element of Eq. (8) to determiner a s . 


In Table I, we compare the FEM value of a s with the experimental value and 

several variational calculations; also included is the Bom-Oppenheimer result. It is 
interesting to note that there remains a discrepancy between the theoretical and 
experimental value of a s , that cannot be accounted for by relativistic corrections. 
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TABLE I. Experimental and theoretical values for the dipole-polarizability of . 


Author (year) 

Ref. 


Bishop and Lamb (1988) Born-Oppenheimer 

[12] 

3.1713 

Jacobson et. al (1997) Experiment 

[8] 

3.168 1(7) 

Shertzer and Greene (1998) 

[13] 

3.168 2(4) 

Bhatia and Drachman (1999) 

[9] 

3.168 0!;™° 

Taylor et al. (1999) 

[14] 

3.168 725 6(1) 

Moss (1999) 

[15] 

3.168 726 

Jacobson et al. (2000) Experiment 

[16] 

3.167 96(15) 

Korobov (2001) 

[17] 

3.168 725 76 

Hilico et al. (2001) 

[18] 

3.168 725 803(1) 

Yan et al. (2003) 

_Q9 L 

3.168 725 802 67(1) 


DIRECT CALCULATION OF THE SCATTERING AMPLITUDE 

We present here a new approach to scattering which does not use partial wave 
analysis [20], The basic idea is to reduce the Schrodinger equation to a set of coupled 
integro-differential equations that can be solved with FEM for the scattering wave 
function. The wave function is then used in the integral expression for the scattering 
amplitude. 

e-H Scattering (static approximation with exchange) 

Our first application is electron-hydrogen scattering in the static approximation 
(with exchange). The trial wave function is given by 

where y/^lr) is an unknown function. We require that the wave function satisfy the 

Schrodinger equation for the two-electron system subject to the asymptotic boundary 
condition 

ikr 

Yt te) e lkrc6sd + // (#)- — (18) 

r 

as r — > oo . Projecting the Schrodinger equation onto <f> u (r 2 ) , we have 


Ate) 


Vf + v 2 — + k 2 -1 

r \ r i n 2 


^tete) =o 


(19) 
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Carrying out the integration over the coordinates of the second electron, we are left with a 
integro-differential equation for the unknown function y/\ : 

[V 2 + F(r) + * 2 V±(F) ±(k 2 + l)fa s (r)(^ s Wt)± Mr)(0u\j^\ri) = O. (20) 

where V(r) = 2e~ 2r (l + }) . 

Although FEM can be applied to integro-differential equations, the last term 
w is problematic because it involves a numerical integration over d 2 r' for 

each value of r . This is computationally prohibitive. To eliminate this problem, we 
introduce a new function 


«,* (?) = (*, ijSiin*) 

(21) 

which satisfies 

V 2 «i(r) = 8 7T</> u (r)y/*(r) 

(22) 

subject to the boundary condition 


a k (?) (fas 1 wt ) 

(23) 


as r — >• oo . Eq. (20) is now replaced by two coupled integro-differential equations 

[v 2 +V(r) + k 2 ]ri(r)± <f> ls ( r )at(r)±(k 2 +l)h s (r)(</> u \rf) = 0 
V 2 a i ± (r)-8^ li (r)^ k ± (?) =0 

which are solved with FEM. The solution yields y/f(r,0), af(r, 6) and //(#) (where 
we have assumed azimuthal symmetry). 

In general, the scattering amplitude obtained directly from the FEM calculation is 
not accurate unless the FEM grid is very large ( r max — > oo) . The results are sensitive to 
the accuracy of the wave function on the boundary, where it is highly oscillatory. 


Integral formula for the scattering amplitude 

In order to reduce the computational effort and improve the accuracy of the 
scattering amplitude, we employ the integral formula for // (0) given by 


/,*(*) = - 7 - k* *.(>;) 

4 TV J 


-+— 


't’UA ,r 2 )d l r x d\ 2 . 


(25) 
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8 


Figure 1 . Results for | // (6) | 2 obtained with FEM are compared with fully converged 
partial wave results. 


Using the ansatz of Eq. (17), we have 


4 n 



2 2 " 

\ e ~ '^ 1 s( r 2 ) 

. r \ r u . 


Vt 0, A ) (t>\s 0 2 ) ± Vl 0 2 , 9 2 ) </>u (>i )]Ar,Ar 2 . (26) 


The accuracy of // (0) now depends on how well the wave function y/ + k (r, 0) is 

represented in the interaction region. In general, the integral expression is not 
particularly useful because it involves a six dimensional integration. However, using our 
definition for a f ( r , 0 ) , we can analytically integrate over four of the six variable to 
obtain 


f;-(O) = -\j o {kr\inO'sme)e- ikr ' co * 0 ' cose [v {r')y/l{r\d')±<l) u {r')a;{r\9')]sin0'r' 1 dr'd9' 


' k 2 + 


- ^<f> u (r) y/f (r' ,0 ') sinO ' r' 2 dr dd ' . 


(27) 
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Figure 2. Results for | f k ( 9 ) | 2 obtained with FEM are compared with fully converged 
partial wave results. 


The final integration is done numerically. 

Using the integral formula, we can obtain extremely accurate and stable results 
with relatively little computational effort. Unlike partial wave analysis, the 
computational effort is independent of energy. In Figs. (1) and (2), we compare our 
results with a fully converged partial wave calculation [21] for the elastic scattering 
amplitude for e-H in the static approximation (with exchange). 

In order to extend this analysis beyond the static exchange approximation, we 
need to include correlation in the wave function, via an explicit dependence on cos 0 n . 
Eventually we plan to include excitation channels in order to obtain the inelastic cross 
sections. 
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